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Abstract 

The recently proposed Einstein molecule approach is extended to compute the free energy of 
molecular solids. This method is a variant of the Einstein crystal method of Frenkel and Ladd[J. 
Chem. Phys. 81, 3188 (1984)]. In order to show its applicability, we have computed the free energy 
of a hard-dumbbells solid, of two recently discovered solid phases of water, namely, ice XIII and 
ice XIV, where the interactions between water molecules are described by the rigid non-polarizable 
TIP4P/2005 model potential, and of several solid phases that are thermodynamically stable for an 
anisotropic patchy model with octahedral symmetry which mimics proteins. Our calculations show 
that both the Einstein crystal method and the Einstein molecule approach yield the same results 
within statistical uncertainty. In addition, we have studied in detail some subtle issues concerning 
the calculation of the free energy of molecular solids. First, for solids with non-cubic symmetry, 
we have studied the effect of the shape of the simulation box on the free energy. Our results 
show that the equilibrium shape of the simulation box must be used to compute the free energy 
in order to avoid the appearance of artificial stress in the system that will result in an increase of 
the free energy. In complex solids, such as the solid phases of water, another difficulty is related 
to the choice of the reference structure. As in some cases there is not an obvious orientation of 
the molecules, it is not clear how to generate the reference structure. Our results will show that, 
as long as the structure is not too far from the equilibrium structure, the calculated free energy 
is invariant to the reference structure used in the free energy calculations. Finally, the strong size 
dependence of the free energy of solids is also studied. 
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I. INTRODUCTION 



Since the pioneering work of Hoover et al.^ determining the free energy of molecular solids 
has been an important area of research .-'^ii^i^'^i^i^ One of the most popular methods to com- 
pute the free energy of solids is the Einstein crystal method, proposed by Frenkel and Ladd 
more than two decades ago.^ In this method, the free energy of a given solid is computed by 
designing an integration path that links the solid to an ideal Einstein crystal with the same 
structure as the real solid, for which the free energy can be analytically computed. This 
method was soon extended to molecular solids.^ In this case, in addition to the springs that 
bound each molecule to its lattice position, springs that keep the particles in the right orienta- 
tion must also be added.-i^ Using this technique, the free energy of several atomic and molec- 
ular solids has been computed.^'^^^^^^^^'i^'i^^^^^^^'^s^'^I^'^^so^^^^^^ 

Quite recently a new method to compute the free energies of solids which was denoted as 
"the Einstein molecule" approach has been proposed.— This method consists of a slight 
modification of the Einstein crystal method. In the Einstein crystal method, the reference 
system is an ideal Einstein crystal with the constraint that the center of mass of the system 
is fixed in order to avoid a quasi-divergence in the integral of the free energy change from 
the real solid to the reference system. This constraint introduces some complexity in the 
method. In particular, the derivation of some terms that contribute to the free energy is 
somewhat involved.™"^ The main idea behind the Einstein molecule approach is that the 
derivation of the analytical expressions can be considerably simplified by fixing the position 
of one molecule instead of fixing the center of mass of the system. The Einstein molecule 
approach has been successfully applied to compute the free energy of the hard-spheres (HS) 
and Lennard- Jones (LJ) face centered cubic (fee) solids. Here it will be shown how it can 
be applied to molecular solids. 

Moreover, even though the Einstein crystal method has been extended to molecular solids 
more than twenty years ago, there are several subtle issues concerning the calculation of the 
free energy that are not clear yet. These difficulties are common to the Einstein crystal and 
Einstein molecule approaches. One of these issues concerns the shape of the simulation box. 
For solids with non-cubic symmetry, prior to the computation of the free energy for a given 
thermodynamic state, the solid structure must be relaxed to obtain the equilibrium unit 
cell corresponding to that thermodynamic state. This is not usually a problem in structures 
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with cubic symmetry, as the equihbrium structure is determined uniquely by the lattice 
parameter a. However, in structures with lower symmetry, it is convenient to first perform a 
simulation in which both the edges and the angles that define the simulation box are allowed 
to relax to the equilibrium structure. This can be achieved by, previously to the free energy 
calculation, performing a Parrinello-Rahman NpT simulation. Other alternative would be 
to perform a simulation at constant volume but where the shape of the simulation box is 
allowed to change, i. e., a variable-shape constant volume (VSNVT) simulation.— '^i'^'^ We 
would like to stress the importance of using the equilibrium structure to compute the free 
energy, otherwise the solid could be under some stress that will lead to an increase of the 
free energy. This has already been noted previously,-!^ but, due to its importance, we believe 
that it is worthy to review this point. 

Another difficulty that one might encounter when computing the free energy of molecular 
solids concerns the reference structure that is used either in the Einstein crystal or in the 
Einstein molecule approaches. In simple solids, in which all the particles exhibit the same 
orientation, this does not pose a problem, as the reference structure is chosen simply as a 
solid where all the particles lie on their lattice positions and are perfectly oriented. However, 
for more complex solids, where not all the molecules exhibit the same orientation, the choice 
of the reference structure might be a subtle issue. This is the case, for example, of some 
solid phases of water that exhibit complex unit cells. In this situation several choices are 
possible. One might choose to build the reference structure by using experimental data to 
obtain the position and orientation of the molecules or, alternatively, one might choose to 
perform an energy minimisation, so that each molecule will be located as to minimise the 
potential energy. Another reasonable choice would be to calculate the average positions and 
orientations at the particular thermodynamic state under study. In view of this ambiguity, 
it is of interest to investigate the effect that one choice or another has on the calculation of 
the free energy. 

Finally, another difficulty arises from the strong size dependence of the free energy of 
solids. In particular, for the fee HS solid, several authors have shown that the free energy 
per particle decreases linearly with 1/N, N being the number of particles in the system.— 
As a consequence, the fluid-solid coexistence point also exhibits a strong size dependence 
(note that the finite size effects on the free energy of the fluid and on the equation of state 
of both phases must also be considered). The size dependence of the fluid-solid coexistence 
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point obtained with the values of the free energy from those works^ is in agreement with 
the coexistence points calculated by Wilding and Bruce^^^ using a completely different 
route, the phase-switch Monte Carlo method.— i^i^ This strong size dependence has also 
been observed for other systems, such as for example, the fee LJ solid.— In this case, the 
situation is more complicated because, in addition to the size dependence of the free energy, 
there is also a dependence on the cutoff of the potential. Both effects must be studied 
separately.— This means that in order to perform a rigorous calculation of the free energy 
of a given solid, the free energy must be computed for different system sizes, so that the 
value of the free energy at the thermodynamic limit can be obtained by extrapolation to 
going to infinity. However, this procedure requires performing many simulations to compute 
the free energy of a solid at just one thermodynamic state. Therefore, it would be useful to 
introduce finite size corrections (FSC), i.e., a simple recipe that would allow one to estimate 
the value of the free energy in the thermodynamic limit from simulations of a system of 
finite size. In a previous paper, we have proposed several FSC whose performances were 
assessed for simple atomic models, namely, the HS and LJ model potentials.— The best 
performance was obtained by the so-called Asymptotic FSC, in which the free energy in the 
thermodynamic limit is estimated from the free energy at a finite size by taking the limit 
when A^ tends to infinity in the expression used to compute the free energy. Depending on 
how this limit was taken, three different variants were proposed, and all of them give quite 
reasonable estimates of the free energy in the thermodynamic limit. However, these results 
might not be general and, therefore, it would be of interest to check whether the FSC work 
well also for molecular solids. 

In this paper we will address all these issues concerning the computation of free energy 
of solids. It is our hope that this will contribute to encourage other authors to compute free 
energies. The paper will be structured in the following way. First, the recently proposed 
Einstein molecule approach will be extended to the case of molecular solids and it will be 
shown that the results obtained for all the solids studied (i.e., a hard-dumbbells solid, a solid 
made of anisotropic particles with octahedral symmetry and the two recently discovered solid 
phases of water, ice XIII and ice XIV) are in agreement, within statistical uncertainty, with 
the results obtained with the Einstein crystal method. Second, the free energy of ices XIII 
and XIV using the rigid non-polarizable TIP4P/2005 model of water will also be calculated. 
These calculations will serve to illustrate the importance of obtaining the equilibrium shape 
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of the simulation box previously to the computation of the free energy and to explore what 
is the best choice for the reference structure that is used in the computation of the free 
energy. Finally, we will perform a systematic study of the size dependence of the free energy 
of several crystalline solids for a simple anisotropic patchy model with octahedral symmetry. 
The performance of the previously proposed FSC will be assessed for this model. 

II. METHOD 

A. Model potentials and solid structures 

In what follows we will consider several pair potentials, for which the intermolecular 
potential will be expressed as: 

N-l N 

Usol = X] UsoliiJ) (1) 

1=1 j=i+l 

where Usoiihj) is the intermolecular potential between molecules i and j. 
1. Hard- dumbbells 

The first model we considered is the hard-dumbbells (HD) model, in which each particle 
consists of two hard-spheres, each of diameter ausi separated by a distance L. The free en- 
ergy of this model has already been studied previously using the Einstein crystal methoct^"^ 
and also theoretically using an extension of the Wertheim theory.— The possible solid 
structures for hard-dumbbells have already been been discussed in previous works.-*^*^ 
Hard-dumbbells can form a hexagonal lattice by arranging the dumbbells in such a way that 
each sphere of a dumbbell lies in a hexagonal layer. The dumbbell axis is then tilted from 
the normal to the layer by an angle equal to arcsin{ ^ ^Vs ^' These layers can be stacked 
as to form a fee lattice (structure designated as CPl) or a hep lattice (structure CP2). In 
these two structures all the dumbbells exhibit the same tilt angle. Another structure can be 
obtained by stacking the layers in such a way that the tilt angle alternates between adjacent 
layers (structure designated as CPS). It has been shown that only the CPl structure is 
thermodynamically stable (for L* = L/cxhs > 0.4).— For values of L* lower than approx- 
imately 0.4, there is a range of pressures for which a plastic fee crystal is the most stable 
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phase.— Finally, it is also possible to form aperiodic fee and hep structures, i.e., structures 
in which the axis of the hard-dumbells are not aligned.— i^i^ The aperiodic fee structure 
becomes thermodynamically stable for values of the elongation L* close to unity.-'^'^'^ As 
the main purpose is to show that both the Einstein molecule approach and the Einstein 
crystal method lead to the same value of the free energy, we have chosen to study only the 
structure designated as CPl for hard-dumbbells with L* = 1 (for this elongation the CPl 
solid is metastable). 

2. The TIP4P/2005 water model. Ices XIII and XIV 

The interaction between water molecules was modelled using a rigid non-polarizable 
model potential, the TIP4P/2005 water model. ^'^ This model is a variant of the TIP4P 
potential,^^ in which the water molecule is modelled by one LJ interaction site on the oxy- 
gen atom, two positive charges located on the hydrogen atom and a negative charge that is 
located on the H-O-H bisector. It has been shown that the TIP4P model is able to predict 
reasonably well the phase diagram of water. It predicts that ice Ih is the most stable solid 
phase at the normal melting point and it reproduces the densities of the solid phases of 
water within 2% of the experimental values.-^ The main failure of this method seems to be a 
melting point about 40K below the experimental value.— It was then clear that the model 
could be improved and several groups proposed variants of this model. In particular, the 
TIP4P/Ice model^ has been fitted to reproduce the experimental melting point of water 
and the TIP4P/Ew^ and TIP4P/2005^ models reproduce the maximum in density at room 
pressure. Among these models, we have chosen to use the TIP4P/2005 model potential, be- 
cause it provides a good description of the phase diagrant^ and also it predicts to good 
accuracy the density of the solid phases of ice.— 

In this work, the free energies of two recently discovered solid phases of water, namely 
ices XIII and XIV,— are computed for the first time. Ice XIII is the proton ordered form 
of ice V. It has a monoclinic unit cell with 28 molecules. Ice XIV is the proton ordered 
form of ice XII. It has a tetragonal unit cell with 12 molecules. The TIP4P/2005 model has 
been shown to reproduce reasonably well the densities of these two solid forms of ice.— In 
the simulations performed in this work the LJ potential was truncated at 8.5 A for both 
solid phases. Standard long range corrections were added to the LJ energy.— Ewald sums 



6 



were used to deal with the long range electrostatic forces. The real part of the electrostatic 
contribution was also truncated at 8.5 A . The screening parameter and the number of 
vectors of reciprocal space considered had to be carefully selected for each crystal phase.— 

3. Model particles with octahedral symmetry 

We also computed the free energy of a patchy model, which has been previously used as 
a simplified model of globular proteins.— "^i^ This model consists of a repulsive core with 
some attractive sites (patches) on its surface. In particular, we studied model particles with 
six patches in an octahedral arrangement. The repulsive core is modelled by the LJ repulsive 
core, while the attractive term is described by the LJ tail modulated by Gaussian functions 
centred at the positions of each patch. Therefore, the total energy between two particles is 
described by the following function: 

Mi J {Tij ) exp ( 1 exp ( - \ r,j > aij 

where ULj{rij) is the Lennard- Jones potential, a is the standard deviation of the Gaussian, 
Ok,ij is the angle formed between patch k (/) on atom i (j) and the interparticle 

vector Yij ( Fj-j), and k^m i^min) is the patch that minimises the magnitude of this angle. 
Additionally, for computational efficiency, the potential is truncated and shifted using a 
cutoff distance of 2.5 au. 

Using reduced units (i.e., choosing the unit of energy and length as the values of the 
LJ parameters Elj and (Tlj), the only parameter that needs to be specified is the width 
of the patches a. In this work, we have chosen a = 0.3 rad., as for this value the whole 
phase diagram has already been studied.— In this previous study, it has been shown that 
there are several solid phases that are thermodynamically stable, namely, simple cubic (sc), 
body-centred cubic (bcc), face-centred cubic (fee) and, at high temperatures, a plastic fee 
crystal. In this work, we will compute the free energy of the three orientationally ordered 
structures (sc, bcc, and fee) for several system sizes. 
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B. The Einstein molecule approach for molecular solids 



As mentioned in the Introduction, the Einstein molecule approach is a variant of the 
Einstein crystal method of Frenkel-Ladd^ that has been proposed quite recently.-22 Analo- 
gously to the Frenkel-Ladd method, the free energy is computed by integration to a reference 
system whose free energy can be computed analytically. The difference is that in the Ein- 
stein molecule approach the reference system is not an ideal Einstein crystal, but an ideal 
Einstein molecule. The Einstein molecule is defined as an Einstein crystal in which one of 
the particles does not vibrate. The name of Einstein molecule has been chosen by analogy 
with molecules, where it is common to use one of the atoms to define the position of a 
molecule, and the vibrational movement of the remaining atoms is given relative to this 
reference atom. The Einstein molecule approach has been successfully applied to compute 
the free energy of simple atomic systems (HS and LJ),'^'^ but we will see that it can be easily 
extended to molecular solids. 

We will start by writing the partition function of a molecular system in the canonical 
ensemble: 

Q = jy|^3jv / (iW[-pU{ri,uJi,...,rN,UN)]driduJi...drNdujN (3) 

where = (xi,yi,Zi) is the position of the reference point of molecule i in Cartesian co- 
ordinates, and uji stands for a set of normalised angles (i.e., J duji = 1) defining the ori- 
entation of particle i. q' = qrqvQe, where qr, qv and qe are the rotational, vibrational, 
and electronic partition functions, respectively. A is the thermal de Broglie wavelength 
(A = (h'^ /(2TTmkBT)Y^'^). There is some freedom in choosing the reference point of the 
molecule. It can be chosen as the center of mass or, alternatively, this reference point can 
be chosen so that all elements of symmetry pass through it (for a more detailed discussion 



see Ref. 



381 ). We have chosen the reference point to be at the center of the sphere for the 



octahedral patchy model and for spheric particles, at the center of mass for hard-dumbbells 
and at the oxygen atom for water. 

The intermolecular potential U depends only on the relative distance between the 
molecules, not on their absolute positions, i.e., it is invariant under translations. This invari- 
ance of the system can be used to write the partition function in a more convenient way by 
performing a change of variables from (ri, r2, rj^j) to (ri, = r2 — ri, r^v' = r^r — ri). 
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Therefore, Equation [3] can be written: 




dri / exp [— /3f/(u;i, r'a, UJ2, r'^, uJn)] diUidr2duj2. ■ .dr'j^duN 



(4) 



The integral k does not depend on the position of particle 1, ri. Therefore, the integration 
over Ti can readily be performed: 



For a system of indistinguishable particles and for a given position of particle 1 there are 
(A^ — 1)! possible permutations of the remaining A^ — 1 particles. The term n can be evaluated 
by computing the integral for a given permutation of the particles [k') and multiplying it 
by the number of permutations, so that the partition function can be written as: 



We will assume that q' has the same value in the two coexisting phases, so that its value 
does not affect the coexistence point. For simplicity, in what follows, we will assign q' the 
value unity. 

We will extend now the definition of the ideal Einstein molecule to molecular solids. For 
atomic solids, an ideal Einstein molecule was defined as an ideal Einstein crystal in which 
one of the particles does not vibrate and acts as reference. For molecular solids, the ideal 
Einstein molecule is defined as an ideal Einstein crystal in which the reference point of 
particle 1 is fixed, but rotations of the molecule about this point are allowed. The reference 
point of particle 1 is called the carrier, because it transports the lattice, i.e., the position of 
the lattice is uniquely defined by the position of the reference point of particle 1. The lattice 
can move as a whole over the volume of the simulation box, and its position is defined by 
the position of the reference point of particle 1. The potential energy of the ideal Einstein 
molecule is given by: 




(6) 





i=2 



i=2 



N 




1=1 
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where is the position of the reference point of molecule i in the reference Einstein solid, 
while Fj represents its position in the current configuration. As can be seen in EqJTJ all 
the particles except particle 1 (which is fixed) are attached to their lattice positions by 
harmonic springs. An orientational field {UEin,or) that forces the particles to adopt the right 
orientation is also included (this field acts over all the particles of the system, including 
particle 1). The orientational field depends on the symmetry of the particles and, thus, an 
orientational field must be defined for each model potential. The orientational field used for 
each one of the model potentials that have been studied in this work will be given in Section 

Kcl 

The partition function of the ideal Einstein molecule can be obtained by performing the 
integral n' for this particular case: 



K 



Ein—mol—id 



exp [-/3A£;(r - ro)^](ir 

3(Af-l)/2 



(7V-1) r 



exp {-(3uEin,or)duJ 



TT 



Q 



Ein,or 



(8) 



where QEin,or is the orientational partition function, which is usually evaluated numerically 
(more details are given the Section Hi CI) . 

The free energy of the ideal Einstein molecule can be obtained by replacing the partition 
function given by Eq. [8] in Eq. (6) 
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(9) 



The numeric value of the thermal de Broglie wavelength A is irrelevant to compute the 
coexistence point as long as the same value is used for both coexisting phases. Therefore, 
we have chosen to assign A the value of the characteristic length for each model potential. 
Thus, for HS A = ansi fo^' HD A = aHS) for LJ A = ctlj, for water A = 1 A and for the 
patchy model A = cr^j. 

In the Einstein molecule approach, the free energy of a given solid is estimated by de- 
signing a path from the ideal Einstein molecule (whose free energy can be computed by Eq. 
in]) to the real solid. This path can be divided into three steps (see Figure [T]). In the first 
step, the ideal Einstein molecule or, what is the same, the position of the reference point 
of the carrier (molecule 1) is constrained to a given position. In the second step, the ideal 
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Einstein molecule with fixed molecule 1 is transformed into the real solid with fixed molecule 
1. Finally, in the last step, the solid of interest is recovered by removing the constraint over 
the position of molecule 1. The free energy change that results from the transformation in 
the first step is given by a term kBTln{V/A^), while the third step contributes by a term 
—kBTln{y/K^). The term V comes from the constraint on the position of molecule 1, and 
the term A'^ comes from the constraint on the momentum. Therefore, the contributions to 
the final free energy of steps one and three cancel out and all what is needed is to compute 
the free energy change between an ideal Einstein molecule and the real solid, both with 
the position (but not the orientation) of particle 1 fixed. This free energy change will be 
computed in two stages. In the first stage we will evaluate the free energy change between 
the ideal Einstein molecule (there is no interaction between the particles, only the external 
Einstein crystal field is present) and the interacting Einstein molecule (in which both the 
springs and the intermolecular potential are present), both with the position of particle 1 
fixed, by a perturbative approach:— 

AAi = Uiattice - fcijT In (exp [-[3{Usol - UiatUce)]) Ein-mol-id ' (l^) 



where Usoi is the potentia. 
frozen lattice (see Ref. 
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energy of the real solid and Uiatuce is the potential energy of the 
for a more detailed discussion). The brackets with the subscript 
Ein — mol — id indicate that the average is performed by sampling the configurations in a 
system where only the Einstein field is present. In the second stage, the interacting Einstein 
molecule with fixed molecule 1 is transformed into the real solid with fixed molecule 1, by 
slowly turning off the springs, according to the following expression: 

U{\) = XUsol + (1 - \){UE^n-n.ol-^d + U ,ol) (H) 

where A is a parameter that takes values between and 1. The free energy change corre- 
sponding to this transformation can be estimated by numerically evaluating the following 
integral: 

= - r ^^"""'r'^ '''''''''' diXAj,). (12) 

Jo 

This integral is usually performed by using a Gauss- Legendre quadrature formula. For that 
purpose, the integrand of this expression must be evaluated at several values of AA^;, which 
can be done by performing NVT MC simulations for those values of the coupling parameter. 
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Taking all the contributions together, the free energy of solid can be computed as: 

-^sol -^Ein—mol—id 

+ AAi + AA2 (13) 

which is the central result of this work. 

An alternative proof of the Einstein molecule approach can be found in the Appendix. 
We show that the Einstein molecule method can be obtained as the limit case of the Einstein 
crystal method when the mass of molecule 1 is much larger than the mass of the remaining 
molecules. 



C. Free energy of the orientational field 



We have said before that the orientational field must be chosen so that it has the same 
symmetry as the molecules. In this section, the orientational fields used for each of the 
studied model potentials are given. In particular, for hard-dumbbells (-Doo,/i symmetry), we 
have chosen the orientational field:- 



N 



UEin,or = ^ [^E,b siu^ {ipb 



(14) 



i=l 



where ipb^i is the angle formed between the axis of particle i and the equilibrium position 
of the axis of particle i in the CPl HD solid. In this case, the partition function of the 
orientational field can be computed as: 



Q 



Ein,or 



1 

An 



exp [—PAEfisin'^{^b,i)) sin 



N 



(15) 



where 6 and (p are the polar angles that define the orientation of the axis of the molecule. In 
this case, the angle can be identified with the polar angle 6. Therefore, this expression 
can be simplified to the following integral in one dimension: 

QEin,or 



exp[l3AE,bix — l)]dx 



(16) 



This integral can be evaluated using a numerical integration method, such as, for example, 
the Simpson's rule. 

The water molecule exhibits symmetry and, therefore, a convenient choice of the 
orientational field is:- 



Ein,or 



N 

E 

1=1 



AE,aSm^ (^a,i) + AeM 



TX 



(17) 
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In this case, the orientation of the molecule is defined by two unitary vectors, a and b. These 
vectors are obtained as the subtraction (a) and the addition (6) of the two hydrogen vectors 
given relative to the position of the oxygen atom. The angle il'a,i is the angle formed by 
the vector a of molecule z in a given configuration {Si) and the vector a in the reference 
configuration (oj q) of the external Einstein field. ipf,^i is defined analogously but with vector 



b (for further details see Ref. l38l ). The orientational partition function can be computed as: 



Q 



Ein, or 



^7r2 



J exp (~/3 < kE,aSin^ 



N 



(18) 



where 6, and x ^ire the Euler angles that define the orientation of the molecule. This 
integral can be simplified by choosing the vector bo as the z axis, so that the Euler angle 
6 is identical to the vector ipf,. It can be evaluated numerically by using a Monte Carlo 
integration method. The efficiency of the Monte Carlo integration method can be consider- 
ably improved by realizing that, for large values of AE,b, the exponential decays very rapidly 
to zero as the angle 6 increases, i.e., as the particle rotates away from the reference orien- 
tation. Therefore, much efficiency is gained by sampling only small values of 6. We have 
chosen to sample cos^ and only those angles for which the cosine is between 0.99 and 1 have 
been considered. About 5000 x 10^ MC cycles were used to evaluate this integral. In a 
previous paper, it has been shown that some approximations can be made to this integral 
for large values of the coupling parameter.- We have found that, for a coupling parame- 
ter AE,a/ {ksT) = AE,b/ {ksT) =25000, the approximation gives a value for the free energy 
of the orientational field, AEin,or/ {NksT) = —l/Nln{QEin.or), about 0.04 lower than that 
obtained by performing the exact integral using the Monte Carlo integral method. In par- 
ticular, using the exact integral we obtained that AEin,or/ {NksT) = 16.05, while using the 
approximate formula, we obtained that AEin,orl {NksT) = 16.01. Although this difference 
is not too large, we recommend to use the exact expression of the integral, using a numerical 
algorithm to evaluate it. 

As with regard to the patchy model with octahedral symmetry (point group Oh), the 
orientational field was:— 

TV 

UEin,or — ^ ^ [A^^^sin {'ipa,i,min 

) + AE,bSm^ {'ipb,i,min)] ■ (19) 

j=l 

where ipa,i,min is the minimum angle formed by any of the vectors that define the position of 
the patches in the particle's reference system with respect to the x axis of a fixed reference 
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system and ipb,i,min is the analogous quantity with respect to the y axis, where the fixed 
reference system has been chosen to be coincident with the orientation of the patches in the 
perfect lattice. Therefore, the orientational partition function is given by: 



Q 



Ein,or 



^ j exp {-P{AE,aSin'^{i'a,min) + ^E,bsin^{i'b,min))} sin9d6d(f)dx 



N 



(20) 



In this case, the integral was evaluated numerically using the Monte Carlo integration 
method and using at least 10^ points. 

In all the cases, we have chosen that both the translational and orientational field have the 
same numeric value of the coupling parameter A^; = A^; ^ = A^; ^. Note, however, that the 
coupling parameter of the translational field, A^; has units of energy over a squared length, 
whereas the orientational coupling parameters A^; „ and A^; {, have dimensions of energy. 

Once the orientational field has been chosen, we can write the explicit form for the integral 
AA2. For example, for water: 



AAn 



N 



N 



EC 



sin^ {ija,i) + ( — 



dA' 



(21) 



N,V,T,A 



i=2 i=l 

where the brackets with the subscript A^, V, T, A' means an average over a simulation of a 
system where both an ideal Einstein field with coupling parameter A' (where A' = AA^;) 
and the solid potential are present (i.e., the total potential is Usoi + A' J2iL2 i^i ~ ^io)"^ + 
a' ^^-^(sin^(?/'(j^j) + (^)^)). For convenience, we will split this expression in two terms, 
one that accounts for the translational contribution {AA2^t) and other that accounts for the 
orientational contribution {AA2^or)'- 
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D. Finite size corrections 



It is well known that the free energy of solids exhibits a strong size dependence.- 
In a recent paper, we have made an attempt to propose some recipes to correct for this 
strong size dependence in a simple way. In what follows, we briefiy review those FSC (a 



more detailed discussion was already given in Ref. 



37|). 
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The two first proposed FSC, namely, the Frenkel-Ladd FSC (FSC-FL) and the hard- 
spheres FSC (FSC-HS) consist of simply adding a term to the free energy of a system of 
particles to obtain an approximation to the free energy in the thermodynamic limit: 

^^sc-FL , , . , Ao..(iV) 2lnN 
NksT ^^^^^>- NksT + N ^^^^ 

NksT NksT +iV ^^^^ 

These are empiric corrections that have been shown to improve the results for the HS fee 
solid. Also we have noted that the term -^InN is approximately equal to the term 7 /N 
except for very small values of N . Therefore, we decided to explore also the performance of 
this FSC: 

AFSC-HS2 , A,oum SlnN 

NksT ^^)- j^ksT + 2 N ^^^^ 

In a second family of FSC which was designated as FSC- Asymptotic, the free energy in 

the thermodynamic limit is estimated by taking the limit when A^ tends to infinity in the 

analytical expression of the free energy of the ideal Einstein molecule (Eq. [H]). Three different 

variants of the FSC-Asymptotic were proposed differing on whether a further approximation 

to the term AA2 was also made. In the first variant {AFsc-asi), no approximation was made 

to AA2: 

Apsc-asi,., , 3 f A^f3AE \ , Ao,,r , AA^{N,Ae) , AA2{N,Ae) 
-^^^(AT^oo)^-/.^— — ^^^^ + ^^^^ (27) 

In a second variant, an approximation to AA2 is made based on the assumption that all 
the A^ — 1 oscillators contribute by the same amount to the integral. This is a reasonable 
approximation for atomic solids (for a fee HS solid with N=108 particles, we obtained 
that all the atoms except the first nearest neighbours contributed approximately by the 
same amount; the contribution of the nearest neighbours is about a 10% lower than the 
contribution of the remaining atoms). For molecular solids, it is important to notice that 
there are two contributions to AA2, one translational and one orient at ional. As pointed out 
before, the Einstein molecule only imposes the constraint on the position of particle 1, but 
not on its orientation. Therefore, assuming that all the molecules contribute by the same 
amount to the translational integral, AA2 can be approximated by the following expression: 

AA2 AA2 t AA2 or N -I ^ AA2 or ( 1 \ . AA2 or , ^ 

' + ■ = It + ' = 1 It + ■ (28) 



NknT NksT NksT N NksT \ Nj ' NksT 
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where AA2^t and AA2^or are given in Eq. [22] and [23l and It is the contribution to the 
translational integral of one single arbitrary particle (under the assumption that all the 
particles contribute by the same amount). We shall assume now that the orientational 
contribution is independent of the system size and that the asymptotic value of AA2^t/NkBT 
is If. In the FSC-as2, AA2 is approximated as: 

. , , AA2,or 

.^N ^00)- It + —— (29) 



Therefore, the FSC-as2 for molecular solids must be slight 
obtained for atomic solids (compare with Eq. 35 of Ref. 



33): 



y modified with respect to that 



AFsc-as2,^, . 3, M2/3As\ Ao,„, AAi(iV,As) 
\N (X)) ^ -In + + 



NksT ' ' 2 \ Ti J NkeT NkeT 

AA2ANAE) ^ AA2,t{N,AE)/{NkBT) 
+ NkBT + (1-1/iV) ^ ^ 

Finally, the last variant is obtained as the mean value of the FSC-asl and FSC-as2: 

ApSC-asS,., , 3, fA^PAE\ , Ao,or , AAi{N,Ae) , AA2,or{N,AE) , 

.(AT ^ 00) ~ -In — + —— + TTT-^^ + TTT-^^ + 



NksT ' ' 2 \ n J NkBT NkBT NkeT 

1 fAA2,t{N, Ae) ^ AA2AN,AE)/iNkBT) \ ^^^^ 



NkBT (1-1/A^) 

Notice that in these expressions AAi and AA2 were obtained by the Einstein molecule 
approach. 



III. RESULTS 



A. The Einstein molecule approach 

Before presenting the results of the free energy calculations with the Einstein molecule 
approach, we will show that fixing one molecule in a solid (in the absence of the Einstein 
field) does not affect the structural properties (due to the translational invariance). For that 
purpose, we computed the radial distribution function in a NVT simulation for an atomic 
system, the HS fee solid, and the site-site radial distribution function for a molecular system, 
the hard-dumbbells CPl solid. We will determine the structure both when one particle is 
fixed and when all the particles are allowed to move. For the HS fee solid, we considered a 
simulation box with =108 particles, so that the possible existence of an inhomogeneity 
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would result in an appreciable change in the radial distribution function. As shown in Fig. [2] 
(a), the radial distribution function is identical regardless of whether one particle is fixed or 
not. As with regard to the hard-dumbbells CPl solid, we considered a simulation box with 
only =32 particles. In this case the center of mass of molecule 1 was fixed but molecule 1 
was allowed to rotate. Our results show that the site-site radial distribution function is again 
identical for a system where all the particles are allowed to move and for a system where 
the center of mass of one of the particles is fixed (see Fig. [2] (b)). Note that it is important 
that the dumbbell with fixed center of mass is allowed to rotate. If molecule 1 is frozen 
at a given orientation, the remaining molecules of the solid will not 'see' all the possible 
orientations of molecule 1. Therefore, all the possible orientations of the fixed molecule 
are not sampled and the fixed particle will introduce an inhomogeneity in the system. We 
checked that this is indeed true by computing also the site-site radial distribution function 
for a system where one particle is not allowed to translate and is not allowed to rotate. In 
this case, it is observed that the value of the site-site radial distribution function at contact 
is affected by the constraint on the orientation of the carrier molecule. In particular, we 
obtained that the value at contact is 5.072 when all the molecules are free to move, which 
is equal (within statistical error) to the value at contact when the position of molecule 1 is 
fixed but it is allowed to rotate, 5.070. However, when molecule 1 is not allowed to translate 
nor to rotate, the contact value of the radial distribution is somewhat lower (5.014), which 
means that the constraint on the orientation introduces an inhomogeneity in the solid. 

The validity of the Einstein molecule approach for molecular solids was checked by com- 
paring the free energies of different molecular solids with those obtained using the Frenkel- 
Ladd Einstein crystal method (as implemented by Poison et al.^'^). At this stage, as the 
purpose was to show that both methods lead to exactly the same results, we performed un- 
usually long simulations in order to reduce the statistical error. In what follows, we describe 
the simulation details for each model. 

For the HD CPl solid, we calculated the free energy for a system with A^ =144 (6x6x4 
unit cells) at a number density p* = pcrjjg = 0.590, where ans is the diameter of each 
one of the hard-spheres of a hard-dumbbell. As the solid is not cubic, we first performed a 
Parrinello-Rahman^i^ NpT MC simulation consisting of 5 x 10^ MC cycles for equilibration 
plus another 5 x 10^ MC cycles for taking averages (a MC cycle is defined as A^ attempts 
to translate or rotate a particle plus one attempt to change the the matrix that defines 
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the simulation box). In agreement with previous calculations, the Parrinello- Rahman NpT 
simulations show that the ratio between the two edges of the unit cell (c/a) is slightly 
different from that at close-packing. Besides the changes in the shape of the simulation box 
it is observed that the orientation of the hard-dumbbells is also different from that at close- 
packing. They change from 9 = 35.26° to ~ 32° and from = 30° to ~ 31°. This has 
already been noted by Vega et al.^ Once we have obtained the equilibrium configuration at 
p* = 0.590, the free energy was calculated by using 16 points to evaluate the integral AA2 by 
the Gauss- Legendre quadrature formula. At each point, the integrand of AA2 was evaluated 
by performing a NVT MC simulation consisting of 2 x 10^ MC cycles for equilibration and 
2 X 10^ MC cycles for taking averages at each value of the coupling parameter. The term 
AAi was calculated in a simulation consisting of 2 x 10^ MC cycles for equilibration and 
5 X 10^ MC cycles for taking averages. The maximum value of the coupling parameter used 
for the free energy calculations was Ke/ [kBT / a'jjg) =4000. 

For the patchy model we considered two system sizes = 125 and = 216. In both 
cases, the free energy was computed at the same thermodynamic state, T* = T / {elj /kB) = 
0.2 and p* = p(j\j = 0.763, where Elj and aij are the parameters of the LJ potential. The 
free energy was evaluated by using 20 points to compute AA2, and at each of those points 
we performed a simulation using 2 x 10^ MC cycles for equilibration and 1 x 10^ cycles for 
taking averages. A maximum value of Ke/ {ksT / a\j) = AE^a/iksT) = AEfi/iksT) =20000 
was used. 

Finally, for ices XIII and XIV (two recently discovered solid phases of water that exhibit 
both oxygen and proton ordering), we computed the free energy at p = 1 bar and T = 
80K. The simulation box contained 3x3x2 unit cells (504 molecules) in the case of ice 

XIII and 3x3x5 unit cells (540 molecules) for ice XIV. As mentioned before, neither 
of these solid phases has cubic symmetry. Ice XIII has a monoclinic unit cell and ice 

XIV has an orthorhombic unit cell. Therefore, previously to the computation of the free 
energy, the solid structure was relaxed to the equilibrium. For ice XIII, we obtained that, 
at p = Ibar and T = 80K, the equilibrium simulation box corresponds to a = 20.39 A, 
b = 22.09 A and c = 28.15 A, and (3 = 109, 6° (a and 7 are equal to 90°). As ice XIV has 
orthorhombic symmetry, only the length of the edges of the box were allowed to fluctuate 
in the simulations, while the angles were kept fixed. In this case, it was obtained that, at 
this thermodynamic state, the length of the edges of the simulation box at equilibrium are 
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a = 24.45 A, 6 = 25.17 A and c = 19.72 A. Once that the equihbrium shape of the simulation 
box was obtained, the positions and orientations of the molecules (i.e, the positions of the 
oxygens and hydrogens) in the reference structure were taken from the crystallographic data 
provided by Salzmann et al.^^ The NpT simulations consisted of 5 x 10^ MC cycles for 
equilibration and 1.5 x 10^ cycles for taking averages. The free energy was computed using a 
maximum value of Ke/ {ksT / A ) = AE,a/ (kBT) = h.Efi/{kBT) =25000. 16 points were used 
to evaluate the term AA2 and, at each of these points, a simulation consisting of 5 x 10^ MC 
cycles (3 x 10^ cycles for equilibration) was carried out, while the term AAi was obtained 
from a simulation consisting of 1 x 10^ MC cycles (2 x 10^ for equilibration). 

The results of the free energies for these systems, as calculated using the Einstein crystal 
method and the Einstein molecule approach are shown in Tables HI and HIl In addition to the 
total free energy, the value of the different terms that contribute to the free energy in both 
methods are also shown. It can be seen that both methods give the same value of the free 
energy within the statistical uncertainty. Although the contribution of the different terms 
that contribute to the free energy is not the same when the center of mass is fixed or when 
molecule 1 is fixed, their sum is invariant. It is interesting to note that the difference between 
AAl/NkeT and AA2/NkBT is about ^InN (see discussion in Ref. l38l ). Therefore, our 
results show that, indeed, the Einstein molecule approach is a valid route to compute the 
free energy of molecular solids. 



B. Free energy of ices XIII and XIV. 

As this is the first time that the free energies of ices XIII and XIV are given, we decided 
to perform more extensive calculations in this case. The effect of the shape of the simulation 
box on the free energy was also studied. Moreover, as mentioned before, it is not obvious 
what orientation of the water molecules should be chosen in the reference structure. For that 
reason, we decided to explore some of the possible orientations to see whether this choice 
affects the results of the free energy calculations. The results presented in this section have 
been obtained from shorter simulations than those of the previous section. Typically each 
simulation consisted of 2 x 10^ MC cycles (4 x 10^ for equilibration). This is usually enough 
to obtain a reasonable accuracy. 

First we calculated the free energy of both phases using the Einstein molecule ap- 
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proach in three different thermodynamic states, namely, p=lbar and T=80K, p=lbar and 
T=250K and p=5000bar and T=80K, so that there are two points along one isobar and 
two points along one isotherm. The values of the free energies at those thermodynamic 
states are shown in Table lllli These data will serve us to check that our calculations 
are thermodynamically consistent, i.e., that the results of the free energy calculations are 
the same as those obtained by thermodynamic integration of the equation of state. For 
ice XIII, through free energy calculations, we obtained Ai{80K,lbar) = — 77. 51(4:) NksT, 
A2{80K,5000bar) = -7 7. 39 (A) NkeT and As{250K,lbar) = -18.51(4:) NksT. Starting 
from Ai and performing thermodynamic integration along this isotherm we obtained that 
^2(807^', 50006ar) = — 77.40(6) NksT. The integration along the isobar starting from Ai 
yields y43(250i^', Ibar) = — 18. 46(6) NkBT. Both estimations are in agreement with the re- 
sults obtained by means of free energy calculations. A good agreement was also obtained for 
ice XIV. In this case, the free energy calculations provide Ai(80K, Ibar) = —77 .82(4) N ksT , 
A2(80K,5000bar) = -77. 73(4) NksT, and A3(250K,lbar) = -18.45(4) NkeT. Starting 
from Ai and integrating along the isotherm 80K, we obtained that A2(80K,5000bar) = 
— 77. 74(6) NksT, and integrating along the isobar Ibar, it is obtained that ^3(2507^, Ibar) = 
— 18. 52(6) NksT, again in good agreement with the results of free energy calculations. 

Once we have confidence on the reliability of our calculations, we studied the effect of the 
shape of the simulation box on the free energy. For that purpose, for ice XIV, the free energy 
was also computed for a simulation box that has been deformed from the equilibrium shape 
at T=80K and p=l bar. The length of the edges at equilibrium L^fl, Ly^ and L^^, were 
deformed to L'^ = L^^ x a, L'y = Ly^ x a, and = Lz^/a"^, so that the density remains 
invariant under this change of the simulation box. Our results show that the free energy 
increases when the shape of the simulation box does not correspond to that at equilibrium. 
In particular, when a deformation defined by a = 0.96 is applied, the free energy increases 
from its value at equilibrium AsoiKNkeT) = -77.82(4) to AsoiKNkeT) = -77.00(4). An 
increase is also found when the edges are scaled with a = 1.04, for which it was found 
that Asoi/ (NksT) = —77.10(4). This is the expected result, as the equilibrium structure 
corresponds to a minimum of free energy, and any perturbation will result in an increase of 
the free energy. These results evidence the importance of obtaining the equilibrium structure 
previously to the computation of the free energy. Otherwise, we would be overestimating 
the value of the free energy. 
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As mentioned before, the positions of the oxygens and hydrogens (i.e., the position and 
orientation of the water molecules) in the reference structure to be used for the Einstein 
field were taken from the experimental values (for ices XIII and XIV both the oxygen and 
hydrogens are ordered). However, the experimental equilibrium positions and orientations 
of the molecules will not be exactly equal to those of the potential model used in the 
simulations. Moreover, it is possible that one would want to study some solid for which there 
are no experimental data available. This does not pose a problem, because, in principle, 
the free energy should not depend on the precise location of the external field provided 
that the field reflects the structure of the system. However, we wanted to check that this 
was indeed true and we computed the free energy using another reference structure. In 
particular, we now used a reference structure in which the water molecules are oriented 
as to minimise the potential energy. This structure was obtained by simulated annealing. 
Starting from a configuration in which the simulation box corresponds to the equilibrium (as 
obtained from the NpT simulation) and where the water molecules have the same positions 
and orientation as those found in experiments, we performed a quenching from 80K to IK, 
using 6 intermediate temperatures, and keeping the shape of the simulation box constant 
along the whole simulation. At each one of the temperatures, the system was allowed to 
evolve during 2 x 10'^ MC cycles. To avoid translations of the system as a whole, we fixed 
the reference point (but not the orientation) of molecule 1 in the annealing. The structure 
obtained from the annealing should be close to the minimum in the potential energy (the 
minimum energy structure would be obtained at OK, at IK the structure is likely to be 
not yet at the minimunt^). Using the structures obtained from simulated annealing, we 
calculated again the free energies of ices XIII and XIV (see Table IIIip . It can be seen 
that the free energy is independent of the reference structure. As expected, the terms AAi 
and AA2 take different values depending on which reference structure has been chosen. 
However, their sum is independent of the reference structure. The term AAi is close to 
the lattice energy and, therefore, it is obvious that its value will depend on the reference 
structure. On the other hand, AA2 is the integral from the real solid to the Einstein molecule 
with intermolecular interactions. In this case, the fact of changing the reference structure 
means that the integral is performed from a new starting point, which results on a different 
value of the integral AA2. However, our results show that the changes in AAi and AA2 
cancel out and, therefore, the same value of the free energy is obtained regardless of which 
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reference structure has been used. Finally we also evaluated the free energy for a reference 
structure in which the atoms are located at their average positions and orientations at that 
thermodynamic state (these were obtained by averaging over 500 snapshots and readjusting 
the positions of the hydrogens to obtain the bond and angle distances of the T1P4P/2005 
model for each molecule), and again the same value of the free energy (within statistical 
error) is obtained. Obviously, it is desirable that the reference structure is close to the 
equilibrium structure, otherwise larger values of the coupling parameter would be needed 
and this will result in a higher statistical error in the evaluation of AA2. 

Taking into account all these considerations, the procedure to compute the free energy 
is schematically described in Fig. [3J It is important that, previously to the computation of 
the free energy, an appropriate reference structure is obtained. 

Before leaving this Section and because this is the first time that the free energy of 
ices XIII and XIV has been computed, we would like to briefly discuss the relative stabil- 
ity of these two ice polymorphs. For that purpose, we computed the chemical potential 
[Pn = {(3A/N) + {P/ pksT)] along the isobars p=lbar and p=5000bar by thermodynamic 
integration (see Fig. H]). It can be seen that, at p =5000bar ice XIV is more stable than ice 
XIII at all temperatures, i.e., from low temperatures up to melting temperatures. On the 
contrary, at p=lbar, ice XIV is slightly more stable than ice XIII at low temperatures, but 
at temperatures close to melting ice XIII seems to be slightly more stable than ice XIV. The 
phase transition seems to occur around T ^ 187K. In any case, it is important to note that 
ice Ih is the most stable phase at p = 1 bar for the TIP4P/2005 model.— 

C. Size dependence of the free energy of molecular solids 

Finally we also studied the size dependence of the free energy of molecular solids by 
analysing the size behaviour of three different solid structures (sc, bcc, and fee) for the 
octahedral patchy model. The free energies of those solid structures as obtained in this 
work using the Einstein molecule approach for several values of are given in Table llVi 
Results for the LJ fee solid at T* = 0.2 and p* = 1.28 are also given in Table IVl In these 
calculations, the LJ potential was truncated at a cutoff distance of 2.7a and long range 
corrections were used (obtained assuming g{r) = 1 beyond the cutoff). Our results show 
that for all the studied solid structures the free energy exhibits a strong size dependence. 
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as found in previous studies.-'^'^i^S'^ It is interesting to note that the slope of the plot 
A{N)/NkBT versus is different depending on the solid structure, even for the same 
model potential. For the patchy model, we obtained that the slope is about -14 for the 
sc structure, about -8 for the bcc and about -12 for the fee. This means that, in order 
to accurately calculate the phase diagram of a given substance, a study of the system size 
dependence must be performed for each considered solid structure, which implies a large 
number of simulations. Therefore it would be useful to have a simple recipe to correct for 
the system size dependence as this could save a large amount of computational time. 

The performance of the FSC was studied for all the considered solid structures of the 
octahedral anisotropic model. In Tables HVl and IVl all the contributions to the free energy are 
explicitly given, as this will allow us to identify the terms that exhibit a stronger dependence 
with the system size. The free energy obtained by applying the proposed FSC to the free 
energy at a given for all the considered solid structures of the octahedral patchy model 
are given in Table IVH and in Figure [H Results of applying the FSC to the calculated free 
energies of the fee L J solid are also given in Table I VII It can be seen that all the proposed 
recipes for finite size corrections give a value of the free energy closer to the thermodynamic 
limit than the estimate obtained from the value of the free energy for a certain A^. The 
FSC-HS, which was based on the slope of the free energy as a function of 1/A^ for HS, 
obviously works better when the slope is similar to the slope of HS (around -7). The same 
is true for the FSC-HS2. Therefore, at a given size, the prediction of the value of the free 
energy in the thermodynamic limit for the sc and fee structures for the patchy model, whose 
slopes were -14 and -12, respectively, is not very accurate. The performance of the FSC-FL 
is also not satisfactory. Although it also seems to give quite good results in some cases (e.g., 
for the bcc solid in the patchy model), there are other solids for which they correct only 
partially for the system size dependence. Finally, the FSC- Asymptotic in its three variants 
seem to give quite accurate results for all the cases studied. This can be understood by 
looking at the size dependence of the terms that contribute to the free energy (see Tables HVl 
andlVj). It can be observed that the terms AAi and the orientational contribution to AA2 
{AA2^or) are almost independent of the system size. All the size dependence comes from the 
free energy of the reference system (as given by Eq. [9]) and, to a much lesser extent, from 
the translational contribution to the integral AA2 {AA2^t)- Therefore, by simply taking the 
limit when A^ ^ 00 in this analytical expression (Eq. [9]), the free energy at a given A^ can 
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be substantially corrected. 

We also calculated the deviation from the correct value (i.e., the free energy in the 
thermodynamic limit) for all the proposed FSC (see Table IVnl) . Data for the HS fee solid 
at three different thermodynamic states taken from a previous work^i are also included in 
Table IVIII The deviation from the correct value {d = ^"'^^ — ^^^^^^zHSi) ^as computed 
for the lowest system size studied in each case. The mean deviation of each FSC computed 
as d = {Yl^=i X 1000 is also given. It can be seen that both the FSC-asl and the 

FSC-as3 exhibit the best performance, obtaining a mean deviation from the correct value 
of 7 or 8 (in 10^^ NksT units). The deviation of the rest of the FSC is not as good, but 
still the mean deviation is typically around 14, which is substantially lower than the mean 
deviation obtained from the true value of the free energy at small values of (around 55). 



IV. CONCLUSIONS 



In this paper, we have extended the Einstein molecule approach to the computation 
of the free energy of molecular solids. The method has been tested using a variety of 
model potentials, which include hard-dumbbells, the TIP4P/2005 water model and a simple 
anisotropic model consisting of a spherical repulsive core with some attractive sites. Our 
results show that both the Einstein crystal method of Frenkel and Ladd (as corrected by 
Poison et al.^^) and the Einstein molecule approach of Vega and Noya give the same results 
of the free energy within statistical accuracy. 

Once the Einstein molecule approach was tested, this method was used to compute the 
free energies of ices XIII and XIV for first time. The free energy was computed at three 
different thermodynamic states, which allowed us to test our free energy calculations by 
performing thermodynamic consistency checks. In addition, we have stressed the importance 
of using the equilibrium shape of the simulation box in the computation of the free energy. 
Our results show that any deformation from this equilibrium structure invariably leads to 
an increase of the free energy. This is the expected behaviour, as the equilibrium structure 
is that that minimises the free energy. Any deformation introduces stress in the system that 
leads to an increase of the free energy. Therefore, for solids with non-cubic symmetry, it is 
important to perform a Parrinello-Rahman NpT simulation to obtain the equilibrium shape 
of the simulation box previously to the computation of the free energy. 
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Moreover, we studied the effect that the choice of the reference Einstein field has on the 
calculation of free energies. In complex solids, such as, for example, ices XIII and XIV, 
there is not an obvious choice of how the water molecules should be oriented in the reference 
structure, as both solids exhibit a complex unit cell with a large number of water molecules 
and in which not all the molecules exhibit the same orientation. We have performed calcu- 
lations of the free energy of both solid phases using a reference structure where the positions 
and orientations of the molecules were taken from experimental data and using a reference 
structure that has been obtained by simulated annealing, i.e., using a reference structure 
that minimises the potential energy (for the equilibrium shape of the simulation box). Our 
results show that, even though the two choices lead to different values of AAi and AA2, the 
addition of both terms is independent of the choice of the reference structure. This is the 
expected result, because we are computing the free energy of the same solid, but using a dif- 
ference reference system (i.e., the position and orientation of the field are slightly different). 
Obviously it is desirable to use a reference structure that is close to the minimum, otherwise 
larger values of the coupling parameter will be needed and this will result in a larger error 
in the evaluation of AA2. This is an important result, because, in many cases, one will be 
interested in real solids with complex structures and the choice of a reference structure will 
be a subtle issue. However, our results show that it is not necessary to obtain the structure 
that minimises the potential energy, as far as the reference structure is not too far from this 
minimum. 

Finally, we have also studied the size dependence of the free energy for a simple anisotropic 
model. Our results show that all the studied solid phases, namely, sc, bcc, and fee, exhibit 
a strong size dependence, although the slope of the plot of A versus 1/N is different for 
each solid phase. In a previous work we also found that the free energy of the fee HS solid 
depends slightly on the density.— This means that there is a complex dependence of the free 
energy with the system size, which depends not only on the model potential, but also on 
the thermodynamic state and on the solid structure. This result seems to suggest that it 
might be difficult to obtain a simple recipe that would allow us to obtain the free energy in 
the thermodynamic limit from the calculated value at a finite size A^. In any case, we tested 
all the previously proposed FSC and we found that the asymptotic FSC-asl and FSC-as3 
manage to give quite accurate estimates of the free energy in the thermodynamic limit for 
all the solids studied so far. 
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Appendix 



We will show that the free energy of the ideal Einstein molecule (Eq. [9]) can be obtained 
as a particular case of the ideal Einstein crystal with fixed center of mass. In the ideal 
Einstein molecule the free energy is given by: 



A 



sol 



Ein—mol~id 



+ AAi + AA2 = Ao + AAi + AA2 



(32) 



The precise expression for Aq is just that of AEm-moi-id (see Eq. [H]). 

In the Einstein crystal method the free energy is computed following the integration path 
shown in Fig. [1] , so that the free energy can be computed as: 



Asol 

A* 



aCM I a a* 



(33) 
(34) 



In this appendix we will show that for a particular choice of the mass of the particles the 
Einstein crystal expression reduces to that of the Einstein molecule expression. The term 
^EiLid is given by : 

A%fn-^d = -kTlniQ%f^,) - kTln{QE^n,ar) (35) 

where QEii,t is given by (see Eq. 97 of Ref. l38|) : 

N 3(W-l)/2 / N \ -3/2 



Mi 



(36) 



■t=i 



and P'"^^(mi, m-Ar) is the contribution of the momenta integral in a system with fixed 
center of mass, where the dependence of P^^{mi, m^r) on the masses is written explicitly. 
The term AA%^ is given by : 



AAI = ksT [ln(P^^^(mi, ...,m;v)/P) - HV/N)] 



(37) 



N 



where -P = 1/(11 ^f) is the contribution to the space of momenta for an unconstrained solid. 

i=l 

Putting together all terms contributing to A^ one obtains: 



Aq = —ksTln 



1 V 



TT 



\j=l 



3{N-l)/2 / N 



-3/2\ 



kBTln{QEi 



(38) 
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Let us now compute Q'e^i the case where all particles of the system have the same 
mass, 777,2 = ^■i = ■■■ = ^N, but the mass of molecule 1 becomes infinitely large compared 
to that of the rest of the particles of the system (this is the choice of Vega and Noya^). In 
this case /ii = 1 and all /i2 = /^s = .. = /iat = (where fii = rrii/ J2iLi "^«)- Obviously under 
these circumstances fixing the center of mass is equivalent to fixing the position of molecule 
1 and, therefore, AA^ = AAi and AA2 = AA2. Let us see if for this particular choice we 
also obtain that Aq = Aq which will complete the proof. Substituting the reduced masses 
/ii = 1 and all /X2 = /^3 = •• = A^iv = in the expression of A^, and after some reordering of 
the terms, one obtains: 

A* = ksTln (^^^ + kBTln ^ ^' - kBTln{QE^n,or) + ksT In (^^^ (39) 

which is exactly equal to the expression of Aq in the Einstein molecule method (Eq. [9]), 
except for the trivial term fc^Tln (^x^^; which obviously will also appear in the fluid phase 
and will not affect the phase equilibria. Thus we have proved that the Einstein molecule 
method can be obtained as a limit case of the Einstein crystal method. 

We have seen that the precise value of P*^*^ is irrelevant to compute the free energy 
(i.e., it does not appear in the expression of A* as given by Eq. [55]) . Nevertheless, for 
completeness, we will compute its value for the two cases we are considering. We will start 
from the general expression of P^^^{mi, ...,171^)'- 



P^*^(777i,..., 777^) = j^^^^ j exp 



N 



i=l 



2mi 



N 



5(2^Pi)cipi...o?Piv (40) 



i=l 



In the particular case that all particles of the system have the same mass (this is the 
choice made by Poison et al.) then for all particles irii = m and = 1/N and then: 

P™(m,...,m) = -^Af-'/^ (41) 



and 

/ ^ x3(iV-l)/2 

See the derivation of this equation in Ref. i38, (Eq. 101). 

Let us now compute P*^^^ for the case where all particles of the system have the same 
mass, rr72 = 7713 = ... = 777 ^r, but the mass of molecule 1 becomes infinitely large compared 
to that of the rest of the particles of the system (this is the choice of Vega and Noya). In 
this case /xi = 1 and all /X2 = /^s = •• = /^at = 0. 
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For this choice of masses: 



1 



N 




P' 



.CM 



/^3(7V-1) 



exp 



5(pi)dpi...rfpjv 



(43) 



i=l 



where Yli=i P* = was simphfied to pi = when the mass of molecule 1 becomes infinitely 
large. 

It is straightforward to integrate this expression to obtain: 
Therefore, the translational contribution to the partition function is: 



which is identical to the expression obtained for the case where all particles have the same 
mass (Eq. HJj). Therefore the expression for A'^fl_^^ is the same when all particles have 
the same mass or for the case where all have the same mass but particle 1 which becomes 
infinitely heavy. 

The partition function of an unconstrained Einstein crystal is given by: 



so that constraining the center of mass in the Einstein crystal amounts to reducing the 
number of degrees of freedom by 3. Notice that this is not the same as A A3 as given by 
Eq. [371 The reason is that Eq. [37] gives the change in free energy for fixing the center of 
mass in a system with translational invariance (i.e., the energy of the system is invariant 
to a translation A of all the particles), and such invariance has been used in the derivation 
leading to the term — ln(V"/A^). Notice that the Einstein crystal does not have translational 
invariance (the energy changes when all the particles are translated by A since the lattice 
does not move), so that AA^ cannot be used to get the free energy change for fixing the 
center of mass in this case. 




(45) 




(46) 
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TABLE I: Free energy of the sc structure of the patchy model particles at T* = 0.2, and of the 
CPl structure of hard dumbbells (HD), as obtained using the Einstein molecule and the Einstein 
crystal methods. For the patchy model we used Ke / {ksT / a\j) = 20000 and A = ctlj and for HD 
AE/iksT/ajjg) = 4000 and A = ans- 





Einstein molecule 




Einstein crystal 


System p* N 


NksT 


NksT 


^sol 

NksT 


A* 

NksT 


AAl AA* A*^, 
NksT NksT NksT 


Patchy (sc) 0.763 125 27.747 


-14.614 


-14.313 


-1.181(7) 


27.689 


-14.614 -14.256 -1.181(7) 


Patchy (sc) 0.763 216 27.792 


-14.614 


-14.311 


-1.134(7) 


27.755 


-14.614 -14.278 -1.138(7) 


HD (CPl) 0.590 144 19.633 


0.001 


-7.056 


12.578(7) 


19.580 


0.001 -7.005 12.576(7) 



TABLE II: Free energies of ices XIII and XIV as calculated using the Einstein crystal and the 
Einstein molecule methods. The simulation box contained A'^ = 504 water molecules for ice XIII 
and N = 540 for ice XIV. Long simulations were performed in order to reduce the statistical error. 
The maximum value of the coupling parameter was 

=25000A-2 and we used A = 1 A. The 
free energy was calculated by performing NVT simulations with the equilibrium simulation box at 
the studied thermodynamic state, namely, T=80K and p=lbar. 







Einstein molecule 




Einstein crystal 


Ice ^?(bar) T(K) p{g/cm^) 


Ao 


AAi AA2 A^^i 


A* 


AAl AAl A,,i 


NksT 


NksT NksT NksT 


NkaT 


NksT NksT NksT 


XIII 1 80 1.262 


29.491 


-91.229 -15.769 -77.508(8) 


29.472 


-91.229 -15.756 -77.512(8) 


XIV 1 80 1.332 


29.493 


-91.073 -16.259 -77.839(8) 


29.475 


-91.073 -16.246 -77.843(8) 
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TABLE III: Free energies of ices XIII and XIV as calculated using the Einstein molecule method. 
The data marked with an asterisk correspond to calculations of the free energy using a reference 
structure in which the positions and orientations of the Einstein field are those obtained from 
simulated annealing up to 1 K, while the data with two asterisks correspond to the structure with 
the average positions and orientations of the water molecules at the particular thermodynamic 
state. As can be seen, the free energy does not depend on the choice of the positions and the 
orientations of the Einstein external field. In all these simulations we have taken A =1 A and 
^E/{kBT/X^) = AE,a/ikBT) = AE,b/ikBT) =25000. 



Ice 


p(bar 


) r(K) p(g/cm3) 


U 






AAi 


AA2 


Asol 


NksT 




NksT 


NksT 


NksT 


NksT 


XIII 


1 


80 


1.262 


-89.08 


25000 


29.49 


-91.23 


-15.77 


-77.51(4) 


XIIP 


1 


80 


1.262 


-89.08 


25000 


29.49 


-92.07 


-14.94 


-77.52(4) 


XIII 


5000 


80 


1.294 


-89.12 


25000 


29.49 


-91.20 


-15.68 


-77.39(4) 


XIII 


1 


250 


1.208 


-26.01 


25000 


29.49 


-28.96 


-19.04 


-18.51(4) 


XIV 


1 


80 


1.332 


-89.64 


25000 


29.49 


-91.07 


-16.24 


-77.82(4) 


XIV* 


1 


80 


1.332 


-89.64 


25000 


29.49 


-92.61 


-14.72 


-77.84(4) 


XIV** 


1 


80 


1.332 


-89.64 


25000 


29.49 


-92.63 


-14.69 


-77.83(4) 


XIV 


5000 


80 


1.360 


-89.71 


25000 


29.49 


-91.02 


-16.20 


-77.73(4) 


XIV 


1 


250 


1.271 


-20.17 


2.5000 


29.19 


-28.95 


-18.99 


-18.-15(4) 
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TABLE IV: Free energies of the patchy model (see Eq. [2]) for different values of N and solid 
structures at T* = 0.2. We also report the value of the three different terms that contribute to 
(see Eq. E]), namely, ^ = ln{A^p)/N, ^ = |/n(A2/3As/^), ^ = -^ln{A^(3AEM. In 
these calculations we used AE/iksT/alj) = 20000 and A was taken as aij- 



System p* N 



Ao,t,l Ao,t,2 Ao,f,3 Ao,„r Aq AAi AAa.t AA2,or A,ai 

NksT NksT NkgT NksT NkgT NkgT NkgT NkgT NkgT 

Patchy (sc) 0.763 125 -0.002 13.138 -0.105 14.716 27.746 -14.614 -5.731 -8.583 -1.181 
Patchy (sc) 0.763 216 -0.001 13.138 -0.061 14.716 27.792 -14.614 -5.728 -8.583 -1.134 
Patchy (sc) 0.763 512 -0.001 13.138 -0.026 14.716 27.828 -14.614 -5.729 -8.582 -1.097 
Patchy (sc) 0.763 1000 -0.000 13.138 -0.013 14.716 27.840 -14.614 -5.729 -8.581 -1.084 
Patchy (bcc) 1.175 250 0.001 13.138 -0.053 14.716 27.802 -13.718 -5.231 -8.562 0.291 
Patchy (bcc) 1.175 432 0.000 13.138 -0.030 14.716 27.824 -13.718 -5.236 -8.564 0.306 
Patchy (bcc) 1.175 1024 0.000 13.138 -0.013 14.716 27.841 -13.718 -5.241 -8.567 0.315 
Patchy (fee) 1.360 256 0.001 13.138 -0.051 14.716 27.804 -6.193 -2.912 -10.108 8.591 
Patchy (fee) 1.360 500 0.001 13.138 -0.026 14.716 27.828 -6.192 -2.913 -10.109 8.614 
Patchy (fee) 1.360 864 0.000 13.138 -0.015 14.716 27.839 -6.190 -2.915 -10.111 8.623 
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TABLE V: Value of the different terms that contribute to the free energy of the LJ fee sohd 
at p* =1.28 and T* =2.0. The LJ potential was truncated at 2.7c7lj. Long range corrections 
(assuming that g(r)=l beyond the cutoff) have been added. We also report the value of the three 
different terms that contribute to j^^'^j. , namely, = l'n{A^p)/N, -^^^ = |/n(A^/3A£;/7r), 

= - J^/n(A2/3A£;/7r). The free energy calculations were performed using a maximum value 
of the coupling parameter KE/{kBT/a\j) = 14000. A was taken as glj- 

\r ^o,t,i ^o,t,2 ^o,t,3 Aq AAi AA2 Asoi 
NksT NksT NksT NksT NksT NksT NksT 



256 0.001 12.603 -0.049 12.555 -3.620 -6.365 2.570 

500 0.000 12.603 -0.025 12.578 -3.620 -6.372 2.586 

864 0.000 12.603 -0.015 12.589 -3.620 -6.377 2.592 

1372 0.000 12.603 -0.009 12.594 -3.620 -6.380 2.594 
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TABLE VI: Free energies of the sc, bcc, and oriented fee crystals for the octahedral patchy model 
at T* = 0.2 and for the LJ model at T* = 2.0 including finite size corrections (FSC). No FSC 
corrections means the true free energy for the system of size N. 



A/{NkT) 



Svqtpm n* N 


No FSC 


FSC-HS2 


FSC-FL 


FSC-HS 


FSC-asl 


FSC-as2 


FSC-as8 




-1.181 


-1.123 


-1.104 


-1.125 


-1 074 


-1.120 


-1 0Q7 


x^atcny l^scj u./oo zio 


-±.1041 


-i.uyo 


1 flK/l 

-l.U04l: 


-l.lUl 


1 n7i 

-l.U ( 1 


1 flQS 

-i.uyo 


-l.UoD 


Patrhv fsrl 763 512 


-1.097 


-1.079 


-1.073 


-1.083 


-1.071 


-1.082 


-1.076 


Patrhv (sc) fl 763 1000 


-1.084 


-1.073 


-1.070 


-1.077 


-1.070 


-1.076 


-1.073 


Patrhv fsrl 768 oo 


-1.069 


-1.069 


-1.069 


-1.069 


-1.069 


-1.069 


-1.069 


jr^dtcny i^Dccj 1.1(0 zou 


n 9Q1 
u.zy 1 


U.OZ4 


U.OOO 




U.o4o 


n "^99 


U.ooz 


Patrhv (hoc) 1 175 482 


0.306 


0.327 


0.334 


0.322 


0.336 


0.324 


0.330 


Patrhv fhrrl 1 1 7'i 1024 

X dtl^llj' \Uyj\^J ±.±tO ±\J^'rt 


81 5 


825 


828 


0.322 


828 


0.322 


825 


Patrhv {hcc^ 1 1 75 no 


0.324 


0.324 


0.324 


0.324 


0.324 


0.324 


0.324 


Patrhv (for") 1 860 2'i6 


8 5Q1 


8 698 


8 684 


8 61 8 


8 641 


8 69Q 


8 685 






o.uoo 


O.UOc/ 








o.uo / 


Patchy (fee) 1.360 864 


8.623 


8.634 


8.638 


8.631 


8.638 


8.634 


8.636 


Patchy (fee) 1.360 oo 


8.637 


8.637 


8.637 


8.637 


8.637 


8.637 


8.637 


LJ 1.28 256 


2.570 


2.602 


2.613 


2.597 


2.618 


2.593 


2.606 


LJ 1.28 500 


2.586 


2.605 


2.611 


2.600 


2.611 


2.598 


2.604 


LJ 1.28 864 


2.592 


2.604 


2.608 


2.600 


2.606 


2.599 


2.603 


LJ 1.28 1372 


2.594 


2.602 


2.605 


2.599 


2.603 


2.598 


2.601 


L.l 1.2<S -x 


2.601 


2.601 


2.601 


2.601 


2.601 


2.601 


2.601 
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TABLE VII: Performance of the different FSC. Deviation of the corrected values of the free energy 
at a given N from the value at the thermodynamic limit (d = — ~^^^^)- For the solids 

studied in this work we computed the mean deviation for the smallest size studied, while for the HS 
solid it was estimated for TV =256. The mean deviation for each FSC is also provided (computed 
as d = Y.\d\/n x lO^). 



System 


P* 


No FSC 


FSC-HS2 


FSC-FL 


FSC-HS 


FSC-Al 


FSC-A2 


FSC-A3 


HS 


1.04086 


-0.028 


0.004 


0.016 


0.000 


0.003 


-0.009 


-0.003 


HS 


1.099975 


-0.030 


0.002 


0.013 


-0.003 


0.007 


-0.008 


0.000 


HS 


1.1500 


-0.034 


-0.002 


0.009 


-0.007 


0.002 


-0.010 


-0.004 
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FIG. 1: Thermodynamic path used in (a) the Einstein molecule approach^ '''^^ and (b) the Einstein 
crystal method.—*^ 
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FIG. 2: (a) Radial distribution function for hard spheres in the fee soUd phases when particle 1 
moves (solid line) and when it does not move (open circles). Results correspond to p* = 1.04086 for 
a system size N = 108. (b) Site-site radial distribution function for hard-dumbbells with L* = 1 
in the CPl structure at p* = 0.590 and for a system size N = 32 when molecule 1 moves (solid 
line) and when the reference point of molecule 1 (i.e., its center of mass) is fixed but molecule 1 
can rotate (open circles). As can be seen, the structural properties are the same, illustrating that 
the properties of the solid present translational invariance. 
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FIG. 3: Schematic representation of the procedure to compute the free energy of molecular solids 
by the Einstein molecule approach. In the first stage a Parrinello-Rahman NpT simulation is 
carried out to obtain the equilibrium shape of the simulation box at the thermodynamic state 
under study. Second, starting from a configuration with the equilibrium shape of the simulation 
box, the position and orientations of the molecules in the lowest energy configuration are obtained 
by simulated annealing (the shape of the simulation box is kept constant during the quenching). In 
order to avoid translations of the system as a whole, the position of the reference point of molecule 
1 is kept fixed during the annealing. The final configuration obtained from this quenching, whose 
energy is Uiatuce, is then used as the reference structure for the computation of the free energy 
(i.e., for the evaluation of terms A^i and AA2)- As described in the text, it is also possible to 
take the structure from the experimental value of the coordinates of the molecules or from the 
average positions. To compute the term A^i, an NVT simulation of the ideal Einstein molecule 
(i.e., with the position of the reference point of molecule 1 fixed) is performed, along which the 
term exp[—P{Usoi — Uiatuce)] is averaged, so that AAi can be computed from Eq. [lOl As with 
regard to the term AA2, several NVT simulations are performed where both the intermolecular 
potential and the Einstein field (for different ■ra.lues of the coupling parameter A') are present, in 
which again molecule 1 is not allowed to translate. For each value of the coupling parameter, the 
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FIG. 4: (a) Chemical potential versus temperature for ices XIII and XIV along the isobars p=lbar 
and p=5000bar at low temperatures, (b) The same as (a) but at temperatures close to the melting 



point of the model. 
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FIG. 5: (Colour online) Size dependence of the free energy of the sc (upper panel), bcc (middle 
panel) and fee (botton panel) solid structures of the six-patches octahedral model. Black circles 
correspond to the true free energy of the system of size N without any FSC correction and the 
black line is a linear fit to these points. The black dashed- lines signals the value of the free energy 
in the thermodynamic limit, obtained from the fit. The red squares correspond to the free energy 
corrected with the FSC-asl, the green diamonds are the free energy corrected with the FSC-as2 
and the blue stars are the values corrected with the FSC-as3. The red, green and blue lines are 
only a guide to the eyes. 
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